function [CONSTANTS IC]=CICR_Constants()
% constants in the model
% global V_m C_m vol_myo vol_JSR  vol_NSR vol_ss A_cap Ca_o R T F CMDN_tot CSQN_tot Km_CMDN Km_CSQN ...
%     v1 v2 v3 tau_tr tau_xfer K_m_up i_CaL_max k_plus_a k_minus_a k_plus_b ...
%     k_minus_b k_plus_c k_minus_c m n E_CaL G_CaL G_CaB K_pcb k_pc_max k_pc_half G_Ca_total...
%      Na_i Na_o k_NaCa K_mNa K_mCa k_sat eta i_pCa_max Km_pCa

V_m = -60;% -82.4202; (millivolt)
C_m = 1; % (microF_per_cm2)
% The dimensions of one regions is 1.5 um (longitudinally) * 3.7 um (transversely) * 2.5
% um (depth)
% The shape of one region is similar to long piece of cake. The length of
% one region is 10 um and the radius is 2.5 um and the angle is 120 (1/3 of
% a full circle)
cakeRadius = 2.5e-6;
cakeLength = 10e-6 ;


volumeConversionFactor = 1e9;
vol_myo = 1/3* cakeRadius^2.0*cakeLength*pi*volumeConversionFactor;%4.5000e-09;%25.84e-6; % (microlitre) -  -  1.5e-6 * 3.7e-6 * 2.5e-6 * (1e9);
vol_JSR =  012 * 1/2 * vol_myo; % (microlitre)   
vol_NSR = 012 * 1/2 * vol_myo; % (microlitre)        
vol_ss = 1.5e-6 * 3.7e-6 * 15e-9 * (1e9); %5.7469e-05*vol_myo;% (microlitre) 
areaConversionFactor = 1e4;
A_cap = 1/3 * 2*cakeRadius*pi* cakeLength *areaConversionFactor; %1.5e-6 * 3.7e-6 * (1e4) ;%1.534e-4; % (cm2) assuming the the demensions are  3.7333 um * 1.5 um - debug
% PMCA constants
i_pCa_max = 2.5;%12;%1;% great for reducing [Ca]i basal level - debug
Km_pCa= 0.5;

% Sodium-related parameters
Na_i  =    14237.1;
Na_o = 140000;
% NCX constants
k_NaCa = 292.8;
K_mNa = 87500;
K_mCa = 1380;
k_sat = 0.1;
eta = 0.35;

Ca_o = 1800; % (micromolar)
R = 8.314; % (joule_per_mole_kelvin)
T = 298;%  (kelvin)
F = 96.5; % (coulomb_per_millimole)

% Calmodulin constants
CMDN_tot = 50; % (micromolar) 
Km_CMDN =0.238; % (micromolar) 

% Calsequestrin   constants
% according to A.Kapla et al paper about smooth
% muscle cells and also the current paper, these numbers(CSQN_tot and Km_CSQN) about Calsequestrin are correct
CSQN_tot =   15000; % (micromolar) 
Km_CSQN =  800 ;% (micromolar) 

% RyR Constants
koCa = 10000;%10; % mM^-2· (ms)^-1; 
kom = 60;%0.06; % (ms)^-1; 
kiCa = 500;%0.5; % (mM)^-1· (ms)-1 ; 
kim = 5;%0.005; % (ms)^-1; 
EC50_SR  = 0.45; %mM; 
ks = 15;%25;%125000;%25%0*10^3;  % (ms)^-1; % NOTE: in Shannon et al the ks = 25 but in Lakatta et al ks = 25e3 ! -  debug
MaxSR = 15; 
MinSR = 1; 
HSR = 2.5; 

% v1 is maximum RyR channel Ca2+ permeability
v1 = 4.5;  % Maximum RyR channel Ca2 permeability (per_millisecond) 
v2 =  1.5*   1.74e-5; %  (per_millisecond)  % it has been changed (modified parameter)
v3 = 0.45;      %  SR Ca2-ATPase (PMCA) maximum pump rate (uM /ms) 

tau_tr = 20;    % (millisecond)
tau_xfer = 8;   % (millisecond)

K_m_up = 0.5; % half-saturation constant for SR Ca ATPase pump 

i_CaL_max = 1;%7;  %(picoA_per_picoF)  - debug
k_plus_a = 0.006075;
k_minus_a = 0.07125;
k_plus_b = 0.00405;
k_minus_b = 0.965;
k_plus_c = 0.009;
k_minus_c = 0.0008;
m = 3;
n = 4;
E_CaL = 63;

% - the problem is that when we add T-type Ca channel, [Ca]i % increases a lot!
G_CaB = 0.01 * 0.000367;%0.28*0.000367; % background Ca conductance   

K_pcb = 0.0005;
k_pc_max = 0.23324;
k_pc_half = 20;

% Initial Condition
P_O1 = 1.498087353057246e-05;%0.149102e-4;     % (dimensionless)
P_O2 =9.596357552441784e-11;% 0.951726e-10;   % (dimensionless)
P_C2 = 1.685348271966624e-04;%0.16774e-3;       % (dimensionless)


O = 6.777812170046545e-12; %0.930308e-18;        % (dimensionless)
C2 = 0.006423125416075;%0.124216e-3;        % (dimensionless)
C3 = 1.557146423337367e-05;%0.578679e-8;        % (dimensionless)
C4 = 1.677756566181471e-08;%0.119816e-12;      % (dimensionless)

I1 = 2.316874514680143e-13;%0.497923e-18;       % (dimensionless)
I2 = 2.716874643061660e-08;%0.345847e-13;       % (dimensionless)
I3 = 7.249477978275548e-08;%0.185106e-13;       % (dimensionless)

Ca_i = 0.115136721404509;     % (micromolar)
Ca_ss = 0.115136721404507;%0.115001;   % (micromolar)
Ca_JSR =  1.302417823198339e+03;%1299.5;     % (micromolar)
Ca_NSR =1.302417827082211e+03;%1299.5;    % (micromolar)


P_RyR = 0; %  (dimensionless)

ryr_R = 0.7499955;
ryr_O = 3.4e-6;
ryr_I = 1.1e-6;
ryr_RI = 0.25;

% ICaL -  L-type Ca channel
G_CaL = (250e-12)*(1e3)/(A_cap*C_m);% (mS/uF)

% ICaT  - T- type Ca channel
G_CaT = (126e-12)*(1e3)/(A_cap*C_m);% (mS/uF)  %0.1729;% 
b_inf = 1.0/(1+ exp(-(V_m + 28.0)/6.1));
g_inf = 1.0/(1+ exp((V_m + 60.0)/6.6));

b = b_inf;
g = g_inf;




CONSTANTS= [ V_m C_m vol_myo vol_JSR  vol_NSR...
    vol_ss A_cap Ca_o R T F CMDN_tot ...
    CSQN_tot Km_CMDN Km_CSQN ...
    v1 v2 v3 tau_tr tau_xfer K_m_up...
    i_CaL_max k_plus_a k_minus_a k_plus_b ...
    k_minus_b k_plus_c k_minus_c m n E_CaL...
     G_CaB K_pcb k_pc_max k_pc_half G_CaT...
     Na_i Na_o k_NaCa K_mNa K_mCa k_sat eta...
     i_pCa_max Km_pCa G_CaL...
     koCa kom kiCa kim EC50_SR ks MaxSR MinSR HSR...
     ];

IC = [O;
    C2;
    C3;
    C4;
    P_C2;
    P_O1;
    P_O2;
    I1;
    I2;
    I3;
    Ca_JSR;
    Ca_NSR;
    Ca_ss;
    Ca_i;
    P_RyR;
    b;
    g;
    ryr_R
    ryr_O
    ryr_I
    ryr_RI
    ];
IC =  IC';
end